results_object = resBAYAREALIKE
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
# States
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.9, splitcex=0.8, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
# Pie chart
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
dev.off()
################ extracting node probabilities
library(dplyr)
#makes list of areas
areas = getareas_from_tipranges_object(tipranges)
areas <- c("AFR", "SA", "NZ", "AUS", "ANT")
numareas = length(areas)
states_list_areaLetters = areas_list_to_states_list_new(areas, maxareas = max_range_size, include_null_range = results_object$inputs$include_null_range)
states=c()
for (i in 1:length(states_list_areaLetters)){
states<-c(states, paste(unlist(states_list_areaLetters[i]),collapse='+'))
}
# marginal probabilities at nodes are here
results_object <- resDEC
probs <- results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
# extracting probabilities of the nodes
nodes<-c(38, 39, 42, 57, 43)
probs = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
probs2 = cbind(nodes,probs[nodes,])
colnames(probs2) <- c("node", states)
probs3 = data.frame(round(probs2, digits = 4))
colnames(probs3)[2] <- "null"
result<-data.frame(matrix(nrow=7, ncol=5))
rownames(result) = c("state1", "state1prob", "state2", "state2prob", "state3", "state3prob", "otherProb")
colnames(result) <- nodes
for (a in 1:5){
x<-as.data.frame(t(probs3[a,2:ncol(probs3)]))
x[,2]<-row.names(x)
colnames(x)<-c("prob", "range")
xx<-arrange(x,desc(prob))
best<-xx[1:3,]
bestTotal<-sum(best[,1])
bestOther<-rbind(best, c(1-bestTotal, "other"))
final = c(bestOther[1,2],bestOther[1,1], bestOther[2,2],bestOther[2,1], bestOther[3,2],bestOther[3,1],
format(round(as.numeric(bestOther[4,1]), digits=4), scientific=FALSE))
result[,a]<-final
}
write.table(result, "nodesPROB_ANTnoCURR_DEC.txt", quote=FALSE, row.names = TRUE, col.names = TRUE, sep="\t")
### DIVALIKE
# marginal probabilities at nodes are here
results_object <- resDIVALIKE
probs <- results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
# extracting probabilities of the nodes
nodes<-c(38, 39, 42, 57, 43)
probs = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
probs2 = cbind(nodes,probs[nodes,])
colnames(probs2) <- c("node", states)
probs3 = data.frame(round(probs2, digits = 4))
colnames(probs3)[2] <- "null"
result<-data.frame(matrix(nrow=7, ncol=5))
rownames(result) = c("state1", "state1prob", "state2", "state2prob", "state3", "state3prob", "otherProb")
colnames(result) <- nodes
for (a in 1:5){
x<-as.data.frame(t(probs3[a,2:ncol(probs3)]))
x[,2]<-row.names(x)
colnames(x)<-c("prob", "range")
xx<-arrange(x,desc(prob))
best<-xx[1:3,]
bestTotal<-sum(best[,1])
bestOther<-rbind(best, c(1-bestTotal, "other"))
final = c(bestOther[1,2],bestOther[1,1], bestOther[2,2],bestOther[2,1], bestOther[3,2],bestOther[3,1],
format(round(as.numeric(bestOther[4,1]), digits=4), scientific=FALSE))
result[,a]<-final
}
write.table(result, "nodesPROB_ANTnoCURR_DIVALIKE.txt", quote=FALSE, row.names = TRUE, col.names = TRUE, sep="\t")
### BAYAREALIKE
# marginal probabilities at nodes are here
results_object <- resBAYAREALIKE
probs <- results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
# extracting probabilities of the nodes
nodes<-c(38, 39, 42, 57, 43)
probs = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
probs2 = cbind(nodes,probs[nodes,])
colnames(probs2) <- c("node", states)
probs3 = data.frame(round(probs2, digits = 4))
colnames(probs3)[2] <- "null"
result<-data.frame(matrix(nrow=7, ncol=5))
rownames(result) = c("state1", "state1prob", "state2", "state2prob", "state3", "state3prob", "otherProb")
colnames(result) <- nodes
for (a in 1:5){
x<-as.data.frame(t(probs3[a,2:ncol(probs3)]))
x[,2]<-row.names(x)
colnames(x)<-c("prob", "range")
xx<-arrange(x,desc(prob))
best<-xx[1:3,]
bestTotal<-sum(best[,1])
bestOther<-rbind(best, c(1-bestTotal, "other"))
final = c(bestOther[1,2],bestOther[1,1], bestOther[2,2],bestOther[2,1], bestOther[3,2],bestOther[3,1],
format(round(as.numeric(bestOther[4,1]), digits=4), scientific=FALSE))
result[,a]<-final
}
write.table(result, "nodesPROB_ANTnoCURR_BAYAREALIKE.txt", quote=FALSE, row.names = TRUE, col.names = TRUE, sep="\t")
}
model = "Antarctica_currents"
x_status="fixed"
x_init = -1.3
x_min = -1.3
x_max = -1.3
x_est = -1.3
areasXXX = c(2,3,5)
for (i in areasXXX){
wd = paste("C:/Users/Martin Fikacek/OneDrive/Manuscripts/MANUSCRIPTS TW/Elmidae VITEK/____new analyses/BioGeoBEARS/optimx_FINAL/X_fixed1.3/", model, "/",i,"/",sep="")  setwd(wd)
for (i in areasXXX){
wd = paste("C:/Users/Martin Fikacek/OneDrive/Manuscripts/MANUSCRIPTS TW/Elmidae VITEK/____new analyses/BioGeoBEARS/optimx_FINAL/X_fixed1.3/", model, "/",i,"/",sep="")
setwd(wd)
extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
extdata_dir
list.files(extdata_dir)
trfn = np(paste(wd, "Elmidae2.newick", sep=""))
tr = read.tree(trfn)
tr
geogfn = np(paste(wd, "Elmidae_geog.txt", sep=""))
tipranges = getranges_from_LagrangePHYLIP(lgdata_fn=geogfn)
tipranges
max(rowSums(dfnums_to_numeric(tipranges@df)))
# Set the maximum number of areas any species may occupy; this cannot be larger
# than the number of areas you set up, but it can be smaller.
max_range_size = i
#######################################################
# Run DEC
#######################################################
# Intitialize a default model (DEC model)
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
# Give BioGeoBEARS the location of the phylogeny Newick file
BioGeoBEARS_run_object$trfn = trfn
# Give BioGeoBEARS the location of the geography text file
BioGeoBEARS_run_object$geogfn = geogfn
# Input the maximum range size
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001    # Min to treat tip as a direct ancestor (no speciation event)
BioGeoBEARS_run_object$include_null_range = TRUE    # set to FALSE for e.g. DEC* model, DEC*+J, etc.
# (For DEC* and other "*" models, please cite: Massana, Kathryn A.; Beaulieu,
#  Jeremy M.; Matzke, Nicholas J.; O'Meara, Brian C. (2015). Non-null Effects of
#  the Null Range in Biogeographic Models: Exploring Parameter Estimation in the
#  DEC Model. bioRxiv,  http://biorxiv.org/content/early/2015/09/16/026914 )
# Also: search script on "include_null_range" for other places to change
# Set up a time-stratified analysis:
# 1. Here, un-comment ONLY the files you want to use.
# 2. Also un-comment "BioGeoBEARS_run_object = section_the_tree(...", below.
# 3. For example files see (a) extdata_dir,
#  or (b) http://phylo.wikidot.com/biogeobears#files
#  and BioGeoBEARS Google Group posts for further hints)
#
# Uncomment files you wish to use in time-stratified analyses:
BioGeoBEARS_run_object$timesfn = "timeperiods.txt"
# BioGeoBEARS_run_object$dispersal_multipliers_fn = "manual_dispersal_multipliers.txt"
#BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt"
#BioGeoBEARS_run_object$areas_adjacency_fn = "areas_adjacency.txt"
BioGeoBEARS_run_object$distsfn = "distances_matrix.txt"
# See notes on the distances model on PhyloWiki's BioGeoBEARS updates page.
# Speed options and multicore processing if desired
BioGeoBEARS_run_object$on_NaN_error = -1e50    # returns very low lnL if parameters produce NaN error (underflow check)
BioGeoBEARS_run_object$speedup = FALSE          # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params)
BioGeoBEARS_run_object$use_optimx = TRUE       # if FALSE, use optim() instead of optimx()
BioGeoBEARS_run_object$num_cores_to_use = 8
# (use more cores to speed it up; this requires
# library(parallel) and/or library(snow). The package "parallel"
# is now default on Macs in R 3.0+, but apparently still
# has to be typed on some Windows machines. Note: apparently
# parallel works on Mac command-line R, but not R.app.
# BioGeoBEARS checks for this and resets to 1
# core with R.app)
# Sparse matrix exponentiation is an option for huge numbers of ranges/states (600+)
# I have experimented with sparse matrix exponentiation in EXPOKIT/rexpokit,
# but the results are imprecise and so I haven't explored it further.
# In a Bayesian analysis, it might work OK, but the ML point estimates are
# not identical.
# Also, I have not implemented all functions to work with force_sparse=TRUE.
# Volunteers are welcome to work on it!!
BioGeoBEARS_run_object$force_sparse = FALSE    # force_sparse=TRUE causes pathology & isn't much faster at this scale
# This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work!
# (It also runs some checks on these inputs for certain errors.)
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
# Divide the tree up by timeperiods/strata (uncomment this for stratified analysis)
BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE)
# The stratified tree is described in this table:
BioGeoBEARS_run_object$master_table
# Good default settings to get ancestral states
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE    # get ancestral states from optim run
# Set up DEC model
# (nothing to do; defaults)
# Look at the BioGeoBEARS_run_object; it's just a list of settings etc.
BioGeoBEARS_run_object
# This contains the model object
BioGeoBEARS_run_object$BioGeoBEARS_model_object
# This table contains the parameters of the model
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","type"]<-x_status
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","init"]<-x_init
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","min"]<-x_min
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","max"]<-x_max
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","est"]<-x_est
# Run this to check inputs. Read the error messages if you get them!
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
res <- bears_optim_run(BioGeoBEARS_run_object)
resDEC <- res
resDEC
resDEC$outputs
resDEC$optim_result$value
calc_loglike_for_optim_stratified
#######################################################
# Run DIVALIKE
#######################################################
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001    # Min to treat tip as a direct ancestor (no speciation event)
BioGeoBEARS_run_object$include_null_range = TRUE    # set to FALSE for e.g. DEC* model, DEC*+J, etc.
# (For DEC* and other "*" models, please cite: Massana, Kathryn A.; Beaulieu,
#  Jeremy M.; Matzke, Nicholas J.; O'Meara, Brian C. (2015). Non-null Effects of
#  the Null Range in Biogeographic Models: Exploring Parameter Estimation in the
#  DEC Model. bioRxiv,  http://biorxiv.org/content/early/2015/09/16/026914 )
# Also: search script on "include_null_range" for other places to change
# Set up a time-stratified analysis:
BioGeoBEARS_run_object$timesfn = "timeperiods.txt"
#BioGeoBEARS_run_object$dispersal_multipliers_fn = "manual_dispersal_multipliers.txt"
#BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt"
#BioGeoBEARS_run_object$areas_adjacency_fn = "areas_adjacency.txt"
BioGeoBEARS_run_object$distsfn = "distances_matrix.txt"
# See notes on the distances model on PhyloWiki's BioGeoBEARS updates page.
# Speed options and multicore processing if desired
BioGeoBEARS_run_object$on_NaN_error = -1e50    # returns very low lnL if parameters produce NaN error (underflow check)
BioGeoBEARS_run_object$speedup = FALSE          # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params)
BioGeoBEARS_run_object$use_optimx = TRUE    # if FALSE, use optim() instead of optimx()
BioGeoBEARS_run_object$num_cores_to_use = 8
BioGeoBEARS_run_object$force_sparse = FALSE    # force_sparse=TRUE causes pathology & isn't much faster at this scale
# This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work!
# (It also runs some checks on these inputs for certain errors.)
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
# Divide the tree up by timeperiods/strata (uncomment this for stratified analysis)
BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE)
# The stratified tree is described in this table:
BioGeoBEARS_run_object$master_table
# Good default settings to get ancestral states
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE    # get ancestral states from optim run
# Set up DIVALIKE model
# Remove subset-sympatry
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","est"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "2-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "ysv*1/2"
# Allow classic, widespread vicariance; all events equiprobable
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","init"] = 0.5
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","est"] = 0.5
# No jump dispersal/founder-event speciation
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = 0.01
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = 0.01
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","type"]<-x_status
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","init"]<-x_init
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","min"]<-x_min
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","max"]<-x_max
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","est"]<-x_est
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
res <- bears_optim_run(BioGeoBEARS_run_object)
resDIVALIKE <- res
#######################################################
# Run BAYAREALIKE
#######################################################
BioGeoBEARS_run_object = define_BioGeoBEARS_run()
BioGeoBEARS_run_object$trfn = trfn
BioGeoBEARS_run_object$geogfn = geogfn
BioGeoBEARS_run_object$max_range_size = max_range_size
BioGeoBEARS_run_object$min_branchlength = 0.000001    # Min to treat tip as a direct ancestor (no speciation event)
BioGeoBEARS_run_object$include_null_range = TRUE    # set to FALSE for e.g. DEC* model, DEC*+J, etc.
# (For DEC* and other "*" models, please cite: Massana, Kathryn A.; Beaulieu,
#  Jeremy M.; Matzke, Nicholas J.; O'Meara, Brian C. (2015). Non-null Effects of
#  the Null Range in Biogeographic Models: Exploring Parameter Estimation in the
#  DEC Model. bioRxiv,  http://biorxiv.org/content/early/2015/09/16/026914 )
# Also: search script on "include_null_range" for other places to change
# Set up a time-stratified analysis:
BioGeoBEARS_run_object$timesfn = "timeperiods.txt"
#BioGeoBEARS_run_object$dispersal_multipliers_fn = "manual_dispersal_multipliers.txt"
#BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt"
#BioGeoBEARS_run_object$areas_adjacency_fn = "areas_adjacency.txt"
BioGeoBEARS_run_object$distsfn = "distances_matrix.txt"
# See notes on the distances model on PhyloWiki's BioGeoBEARS updates page.
# Speed options and multicore processing if desired
BioGeoBEARS_run_object$on_NaN_error = -1e50    # returns very low lnL if parameters produce NaN error (underflow check)
BioGeoBEARS_run_object$speedup = FALSE          # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params)
BioGeoBEARS_run_object$use_optimx = TRUE    # if FALSE, use optim() instead of optimx()
BioGeoBEARS_run_object$num_cores_to_use = 8
BioGeoBEARS_run_object$force_sparse = FALSE    # force_sparse=TRUE causes pathology & isn't much faster at this scale
# This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work!
# (It also runs some checks on these inputs for certain errors.)
BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)
# Divide the tree up by timeperiods/strata (uncomment this for stratified analysis)
BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE)
# The stratified tree is described in this table:
BioGeoBEARS_run_object$master_table
# Good default settings to get ancestral states
BioGeoBEARS_run_object$return_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE
BioGeoBEARS_run_object$calc_ancprobs = TRUE    # get ancestral states from optim run
# Set up BAYAREALIKE model
# No subset sympatry
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","est"] = 0.0
# No vicariance
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","est"] = 0.0
# No jump dispersal/founder-event speciation
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free"
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = 0.01
# BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = 0.01
# Adjust linkage between parameters
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "1-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/1"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "1-j"
# Only sympatric/range-copying (y) events allowed, and with
# exact copying (both descendants always the same size as the ancestor)
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","init"] = 0.9999
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","est"] = 0.9999
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","type"]<-x_status
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","init"]<-x_init
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","min"]<-x_min
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","max"]<-x_max
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["x","est"]<-x_est
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table
# Check the inputs
check_BioGeoBEARS_run(BioGeoBEARS_run_object)
res <- bears_optim_run(BioGeoBEARS_run_object)
resBAYAREALIKE <- res
#########################################################################
#########################################################################
#########################################################################
#
# CALCULATE SUMMARY STATISTICS TO COMPARE
# DEC, DEC+J, DIVALIKE, DIVALIKE+J, BAYAREALIKE, BAYAREALIKE+J
#
#########################################################################
#########################################################################
#########################################################################
#########################################################################
#########################################################################
#########################################################################
# REQUIRED READING:
#
# Practical advice / notes / basic principles on statistical model
#    comparison in general, and in BioGeoBEARS:
# http://phylo.wikidot.com/advice-on-statistical-model-comparison-in-biogeobears
#########################################################################
#########################################################################
# Set up empty tables to hold the statistical results
restable = NULL
teststable = NULL
#######################################################
# Statistics
#######################################################
res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resDEC, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resDIVALIKE, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
res3 = extract_params_from_BioGeoBEARS_results_object(results_object=resBAYAREALIKE, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
resALL <- rbind(res1, res2, res3)
restable <- cbind(model=c("DEC", "DIVALIKE", "BAYAREALIKE"), resALL)
# With AICs:
AICtable = calc_AIC_column(LnL_vals=restable$LnL, nparam_vals=restable$numparams)
restable = cbind(restable, AICtable)
restable_AIC_rellike = AkaikeWeights_on_summary_table(restable=restable, colname_to_use="AIC")
restable_AIC_rellike = put_jcol_after_ecol(restable_AIC_rellike)
restable_AIC_rellike
# With AICcs -- factors in sample size
restable2 = restable
samplesize = length(tr$tip.label)
AICtable = calc_AICc_column(LnL_vals=restable$LnL, nparam_vals=restable$numparams, samplesize=samplesize)
restable2 = cbind(restable2, AICtable)
restable_AICc_rellike = AkaikeWeights_on_summary_table(restable=restable2, colname_to_use="AICc")
restable_AICc_rellike = put_jcol_after_ecol(restable_AICc_rellike)
restable_AICc_rellike
# Also save to text files
write.table(restable_AIC_rellike, file="restable_AIC_rellike.txt", quote=FALSE, sep="\t")
write.table(restable_AICc_rellike, file="restable_AICc_rellike.txt", quote=FALSE, sep="\t")
# Save with nice conditional formatting
write.table(conditional_format_table(restable_AIC_rellike), file="restable_AIC_rellike_formatted.txt", quote=FALSE, sep="\t")
write.table(conditional_format_table(restable_AICc_rellike), file="restable_AICc_rellike_formatted.txt", quote=FALSE, sep="\t")
#######################################################
# Plot ancestral states - make PDF
#######################################################
pdf("ancestral_reconstructions.pdf", width = 20, height = 10)
par(mfrow=c(1,2))
analysis_titletxt ="DEC"
results_object = resDEC
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
# States
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.9, splitcex=0.8, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
# Pie chart
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
analysis_titletxt ="DIVALIKE"
results_object = resDIVALIKE
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
# States
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.9, splitcex=0.8, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
# Pie chart
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
analysis_titletxt ="BAYAREALIKE"
results_object = resBAYAREALIKE
scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
# States
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.9, splitcex=0.8, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
# Pie chart
plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
dev.off()
################ extracting node probabilities
library(dplyr)
#makes list of areas
areas = getareas_from_tipranges_object(tipranges)
areas <- c("AFR", "SA", "NZ", "AUS", "ANT")
numareas = length(areas)
states_list_areaLetters = areas_list_to_states_list_new(areas, maxareas = max_range_size, include_null_range = results_object$inputs$include_null_range)
states=c()
for (i in 1:length(states_list_areaLetters)){
states<-c(states, paste(unlist(states_list_areaLetters[i]),collapse='+'))
}
# marginal probabilities at nodes are here
results_object <- resDEC
probs <- results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
# extracting probabilities of the nodes
nodes<-c(38, 39, 42, 57, 43)
probs = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
probs2 = cbind(nodes,probs[nodes,])
colnames(probs2) <- c("node", states)
probs3 = data.frame(round(probs2, digits = 4))
colnames(probs3)[2] <- "null"
result<-data.frame(matrix(nrow=7, ncol=5))
rownames(result) = c("state1", "state1prob", "state2", "state2prob", "state3", "state3prob", "otherProb")
colnames(result) <- nodes
for (a in 1:5){
x<-as.data.frame(t(probs3[a,2:ncol(probs3)]))
x[,2]<-row.names(x)
colnames(x)<-c("prob", "range")
xx<-arrange(x,desc(prob))
best<-xx[1:3,]
bestTotal<-sum(best[,1])
bestOther<-rbind(best, c(1-bestTotal, "other"))
final = c(bestOther[1,2],bestOther[1,1], bestOther[2,2],bestOther[2,1], bestOther[3,2],bestOther[3,1],
format(round(as.numeric(bestOther[4,1]), digits=4), scientific=FALSE))
result[,a]<-final
}
write.table(result, "nodesPROB_ANTnoCURR_DEC.txt", quote=FALSE, row.names = TRUE, col.names = TRUE, sep="\t")
### DIVALIKE
# marginal probabilities at nodes are here
results_object <- resDIVALIKE
probs <- results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
# extracting probabilities of the nodes
nodes<-c(38, 39, 42, 57, 43)
probs = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
probs2 = cbind(nodes,probs[nodes,])
colnames(probs2) <- c("node", states)
probs3 = data.frame(round(probs2, digits = 4))
colnames(probs3)[2] <- "null"
result<-data.frame(matrix(nrow=7, ncol=5))
rownames(result) = c("state1", "state1prob", "state2", "state2prob", "state3", "state3prob", "otherProb")
colnames(result) <- nodes
for (a in 1:5){
x<-as.data.frame(t(probs3[a,2:ncol(probs3)]))
x[,2]<-row.names(x)
colnames(x)<-c("prob", "range")
xx<-arrange(x,desc(prob))
best<-xx[1:3,]
bestTotal<-sum(best[,1])
bestOther<-rbind(best, c(1-bestTotal, "other"))
final = c(bestOther[1,2],bestOther[1,1], bestOther[2,2],bestOther[2,1], bestOther[3,2],bestOther[3,1],
format(round(as.numeric(bestOther[4,1]), digits=4), scientific=FALSE))
result[,a]<-final
}
write.table(result, "nodesPROB_ANTnoCURR_DIVALIKE.txt", quote=FALSE, row.names = TRUE, col.names = TRUE, sep="\t")
### BAYAREALIKE
# marginal probabilities at nodes are here
results_object <- resBAYAREALIKE
probs <- results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
# extracting probabilities of the nodes
nodes<-c(38, 39, 42, 57, 43)
probs = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
probs2 = cbind(nodes,probs[nodes,])
colnames(probs2) <- c("node", states)
probs3 = data.frame(round(probs2, digits = 4))
colnames(probs3)[2] <- "null"
result<-data.frame(matrix(nrow=7, ncol=5))
rownames(result) = c("state1", "state1prob", "state2", "state2prob", "state3", "state3prob", "otherProb")
colnames(result) <- nodes
for (a in 1:5){
x<-as.data.frame(t(probs3[a,2:ncol(probs3)]))
x[,2]<-row.names(x)
colnames(x)<-c("prob", "range")
xx<-arrange(x,desc(prob))
best<-xx[1:3,]
bestTotal<-sum(best[,1])
bestOther<-rbind(best, c(1-bestTotal, "other"))
final = c(bestOther[1,2],bestOther[1,1], bestOther[2,2],bestOther[2,1], bestOther[3,2],bestOther[3,1],
format(round(as.numeric(bestOther[4,1]), digits=4), scientific=FALSE))
result[,a]<-final
}
write.table(result, "nodesPROB_ANTnoCURR_BAYAREALIKE.txt", quote=FALSE, row.names = TRUE, col.names = TRUE, sep="\t")
}
